GEOPHYSICAL RESEARCH LETTERS, VOL. ???, XXXX, DOL10.1029/, 



GPS radio occultation with GRACE: Atmospheric 
profiling utilizing the zero difference technique 

G. Beyerle, T. Schmidt, G. Michalak, S. Heise, J. Wickert, and Ch. Reigber 
GeoForschungsZentrum (GFZ) Potsdam, Germany 



G. Beyerle, S. Heise, G. Michalak, Ch. Reigber, T. Schmidt, J. Wickert, GeoForschungsZen- 
trum Potsdam (GFZ), Department 1, Geodesy and Remote Sensing, Telegrafenberg, D- 
14473 Potsdam, Germany (e-mail: gbeyerle@gfz-potsdam.de) 



DRAFT 



February 2, 2008, 4:32am 



DRAFT 



X - 2 BEYERLE ET AL.: GPS OCCULTATION USING ZERO DIFFERENCING 

Abstract. Radio occultation events recorded on 28-29 July 2004 by a GPS 
receiver aboard the GRACE-B satellite are analyzed. The stability of the re- 
ceiver clock allows for the derivation of excess phase profiles using a zero dif- 
ference technique, rendering the calibration procedure with concurrent ob- 
servations of a reference GPS satellite obsolete. 101 refractivity profiles ob- 
tained by zero differencing and 96 profiles calculated with an improved sin- 
gle difference method are compared with co-located ECMWF meteorolog- 
ical analyses. Good agreement is found at altitudes between 5 and 30 km with 
an average fractional refractivity deviation below 1% and a standard devi- 
ation of 2-3%. Results from end-to-end simulations are consistent with these 
observations. 
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1. Introduction 

In recent years atmospheric sounding by space-based Global Positioning System (GPS) 
radio occupation (RO) is considered a valuable data source for numerical weather pre- 
diction and climate change studies. From 1995 to 1997 the proof-of-concept GPS/MET 
mission performed a series of successful measurement campaigns [Rocken et al, 1997]; 
since 2001 a RO instrument operates aboard the CHAMP (CHAllenging Minisatellite 
Payload) [Reigber et al, 2004] satellite. 

During an occultation event the RO receiver aboard the low Earth orbiting (LEO) 
spacecraft records dual-frequency signals transmitted from a GPS satellite setting behind 
the horizon. Voltage signal-to-noise ratios (SNR V ) and carrier phases at the two GPS 
L-band frequencies f\ = 1.57542 GHz (LI) and f 2 = 1.2276 GHz (L2) are tracked with 
a sampling frequency of typically 50 Hz. The ionosphere and neutral atmosphere induce 
optical path length deviations from the geometrical distance between transmitter and 
receiver. These deviations are denoted as excess phase paths x an d are related to the 
ray bending angle a as a function of impact parameter p. From a{p) the atmospheric 
refractivity profile N(z) = (n(z) — 1) • 10 6 is derived, where n(z) denotes the real part of 
the atmospheric refractive index at altitude z [Kursinski et al, 1997; Hajj et al, 2002]. 

GPS/MET observations were analyzed with a double difference method to correct for 
clock errors of the GPS transmitters and the receiver aboard the LEO satellite [Rocken 
et al, 1997]. Since deactivation of Selective Availability (S/A), an intentional degradation 
of broadcast ephemeris data, on 2 May 2000, the GPS clock errors are reduced by orders of 
magnitude. Without S/A GPS clocks are sufficiently stable to replace double differencing 
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by the single difference technique thereby eliminating the need for concurrent high-rate 
ground station observations [Wickert et al, 2002]. 

On 28-29 July 2004 the GPS RO receiver aboard the GRACE-B satellite, the second 
of the two GRACE satellites [Tapley et al, 2004], was activated for 25 hours to test the 
receiver's occultation mode [Wickert et al, 2005]. The GRACE-B receiver clock time 
was adjusted to GPS coordinate time only once during the measurement period. Thus, 
the GRACE-B clock solutions can be interpolated from the precise orbit determination's 
(POD) temporal resolution of 30 s to the occultation receiver's resolution of 20 ms and 
the excess phase path is obtained without simultaneous observations of a reference GPS 
satellite (zero differencing). 

2. Methodology and Data Analysis 

The LI and L2 signals transmitted by the occulting GPS satellite are recorded by 
the receiver's signal tracking phase-locked loops. Their numerically-controlled oscillators 
(NCOs) are phase-locked to the signal carriers transmitted by the GPS satellite which in 
turn are phase-locked to the GPS transmitter clock. We may regard the NCOs as clocks 
that are synchronized to the transmitter clock of the occulting satellite. The occulting 
GPS clock time, as recorded by the receiver NCO, is denoted by i°. (k — 1 for LI, k — 2 
for L2; though, in the following the subscript k is omitted.) During an occultation event 
the receiver measures both i° and the corresponding LEO clock time i L . 

2.1. Zero and Single Differencing Methods 

Since the GPS and LEO satellites move within Earth's gravitational potential, the 
proper time t, the coordinate time t and the clock time t in general will not agree. In 
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order to extract the excess phase paths x(t) from the measurements L OL = c{t L — t°), 
the LEO clock corrections (i L — t L ) and the occulting GPS clock corrections (i° — t°) 
have to be taken into account [Hajj et ai, 2002]. Here, c denotes the velocity of light in 
vacuum. We note that both, (i L — t L ) and (t° — t°), include relativistic effects relating 
coordinate time to proper time. Thus, the zero differencing (ZD) equation is 

t L -t° = - l ol - (t L - t L ) - (t° - i°) . (i) 

c 

The left hand side of Eqn. 1 contains the excess phase path x through 

t L - t° « \ (d OL + x) (2) 

where contributions from the gravitational time delay are neglected [Hajj et ai, 2002] and 
d OL = \r° — r L \. The positions of the occulting GPS and the LEO satellite, r° and f L , 
and their clock corrections, (i° — t°) and (i L — t L ), are obtained from the POD [Konig 
et ai, 2004]. CHAMP's and GRACE-B's orbits and clock solutions are provided at a 
temporal resolution of 30 s. 

GRACE-B's receiver clock is adjusted to GPS coordinate time only once during the 
25 hour measurement period on 28 July 16:24 UTC. The difference between clock time 
and coordinate time, (i L —t L ), at the required temporal resolution of 20 ms could therefore 
successfully be interpolated from the 30 s clock solutions. The receiver clock aboard 
CHAMP, on the other hand, is adjusted once per second to achieve a 1 /is maximum 
deviation with respect to coordinate time. (See discussion of Fig. 2 below.) These 
adjustments introduce discontinuities that cannot be corrected for on the basis of the 30 s 
clock solutions and require the simultaneous observation of a referencing GPS satellite 
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(single differencing (SD)) 

t L - t° = - {L OL - L RL ) + {t L - t R ) + (t R - t R ) - (t° - t°) (3) 

with superscript R indicating the reference satellite [Wickert et al, 2002]. The reference 
link observation corrected for ionospheric signal propagation, L RL , is calculated from 

L RL = Cl L%-c 2 L%. (4) 

with parameters c\ and c 2 defined as q = (/i) 2 /[(/i) 2 — (/2) 2 ], i — 1,2. 

GFZ's operational RO processing system [Wickert et al, 2004] improves the ionospheric 
correction of the reference link observations (Eqn. 4) by assuming that quiet ionospheric 
conditions prevail (modified single differencing, in the following abbreviated as SDM). 
Within the SDM approach L RL is replaced by 

L RL = L R [ + c 2 (L*f - L^> (5) 

where (•) denotes a smoothing (low-pass) filter. Since the noise level of L R % is much larger 
than the noise level of L R f (see Fig. 1) the filter effectively reduces the L2 noise. 

2.2. Comparison with ECMWF 

The observed GRACE refractivity profiles are intercompared with meteorological anal- 
ysis results provided by the European Centre for Medium-Range Weather Forecasts 
(ECMWF). Pressure and temperature values from the ECMWF fields are calculated by 
linear interpolation between grid points (0.5° x 0.5° resolution). Linear interpolation in 
time is performed between 6 h analyses fields. The comparison between RO observation 
and meteorological analysis is performed on the 60 ECMWF model levels ranging from 
the ground surface up to 0.1 hPa (about 60 km altitude) after converting the pressure lev- 
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els to height coordinates. Vertical spacing of the model grid points increases from about 
200 m at 1 km altitude to about 700 m at 10 km altitude. 

2.3. Simulation Study 

We interpret the GRACE-B observations using results from end-to-end simulations. 
The refractivity field is assumed to be spherically symmetric and calculated from the 
corresponding ECWMF profile; electron densities are modelled using the Parameterized 
Ionospheric Model (PIM) [Daniell et al, 1995]. The occultation events are simulated with 
circular approximations of the true GRACE-B and GPS satellite orbits. Signal amplitudes 
and carrier phases on the occultation link are calculated from bending angles with the 
inverse CT technique [Gorbunov, 2003] modified for circular orbits (inverse FSI or CT2 
method). The amplitude is normalized to the observation averaged over the first 10 s of 
the occultation. The reference link data are assumed to be multipath-free and the carrier 
phases are determined using geometric optics, their signal amplitudes are taken to be 
constant. Again, their magnitude is extracted from the observation. On both links ray 
bending due to the ionosphere and neutral atmosphere is taken into account. 

In principle, carrier tracking loop errors as well as clock instabilities and ionospheric 
disturbances contribute to the uncertainties of the observations L OL and L RL . For a given 
amplitude SNR V the tracking loop phase noise a c can be estimated from [Ward, 1996] 



with carrier wavelength A, sampling time T and carrier loop bandwidth B w . 

The simulated LI phase noise is calculated from Eqn. 6 using T = 20 ms and B w = 
20 Hz, since this parametrization is in approximate agreement with the observed Lff 
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noise levels which are plotted in Fig. 1 (circles). The full line marks the result from 
Eqn. 6. The corresponding L2 noise values, however, estimated from Eqn. 6 significantly 
underestimate the observations since additional losses from semi-codeless or codeless L2 
tracking are not taken into account by Eqn. 6. A heuristic parametrization is therefore 
used to simulate the L2 phase noise (dashed line in Fig. 1). 

The uncertainties of both observations, Lff and Lf^ , exhibit a distinct dependency on 
SNR V (Fig. 1). Since the clock error is expected to be independent of SNR V we conclude 
from Fig. 1 that tracking loop errors dominate the clock uncertainties. Enhanced noise 
values in some of the observations are presumably attributed to ionospheric disturbances. 
From the simulated SNR v (t) and x(t) data bending angle profiles are reconstructed using 
the FSI method [Jensen et al, 2003] and (above 20 km altitude) using the geometric 
optics method. Finally, refractivity profiles are retrieved by Abel-transforming [Fjeldbo 
et al, 1971] the bending angle profiles thereby closing the simulation loop. 

3. Discussion 

Between 6:03 UTC on 28 July 2004 and 7:09 UTC on 29 July 2004 the GPS receiver 
aboard GRACE-B was activated to test occultation measurement mode. During these 
25 hours the receiver recorded 109 occultation events lasting longer than 40 s; 101 events 
could be successfully converted to atmospheric refractivity profiles with ZD and SDM 
processing, 8 were removed as outliers. A profile is regarded as outlier if the fractional 
refractivity error with respect to ECMWF exceeds 10% at altitudes above 10 km. The 
corresponding yield for SD is 96 profiles. 
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From the reference link observations L RL the difference between receiver clock time and 
proper time, (t L — i L ), can be estimated assuming quiet ionospheric conditions [Wickert 
et al, 2002]. Its temporal evolution, A(t L -t L )/T, extracted from the first GRACE-B oc- 
cultation measurement on 28 July 2004, 6:10 UTC at 55.3°N, 22.3°E, is plotted in Fig. 2, 
bottom panel, black line. Here, A(x n ) = x n+ i—x n denotes the forward difference operator 
and the sampling rate is T = 20 ms. The mean clock drift of 31.3347±0.6305 ns/s is con- 
sistent with a 31.3338 ns/s mean drift obtained from GRACE-B precise orbit calculations 
(white line). 

For comparison the top panel of Fig. 2 shows an occultation event recorded by CHAMP 
on the same day, 0:18 UTC at 69.4°N, 12.5°W. The CHAMP clock drift (-1.0642 ± 
8.5187 ns/s) exhibits periodic discontinuities of about 10 ns/s about every 18 s; in addition, 
the CHAMP clock is adjusted every second in order to meet the design specification of 
1 /is absolute time signal error with respect to coordinate time. These discontinuities 
cannot be modelled from the 30 s POD data eliminating at present the possibility of zero 
differencing with CHAMP. Since the data sets shown in Fig. 2 are extracted from reference 
link observations, the random components contain contributions from the ionosphere as 
well as the tracking loops and should not be interpreted as clock noise. 

The mean fractional refractivity error {(N ohs — N met )/N met ) derived from the GRACE-B 
observations is plotted in Fig. 3 (left panel), where N t, s and N met denote the observed 
and ECMWF refractivities, respectively. The 1-a standard deviations are marked as thin 
lines. The number of extracted data points as a function of altitude is plotted in the right 
panel. 
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The ZD results are characterized by smaller standard deviations at altitudes below 
20 km since noise contributions from the reference link are present in the SD, but not in 
the ZD excess phases. Above 20 km multipath is assumed to be negligible and bending 
angles are calculated using geometric optics. In addition, the SDM results are marked by 
crosses in the left panel of Fig. 3. The almost perfect agreement with the ZD profiles, both 
in terms of bias (crosses) and standard deviation (not shown), emphasizes the impact of 
the respective ionospheric correction procedure (Eqn. 4 vs. Eqn. 5) on the tropospheric 
retrieval results. 

Fig. 3 suggests that ZD processing reduces the lower tropospheric refractivity bias 
compared to the SD results. At low SNR V substantial phase noise contributions are present 
in the excess phase path profile x(t) thereby degrading the radioholographic interference 
patterns encoded in x(0- These patterns originate from interfering rays associated with 
large negative vertical refractivity gradients that are caused by humidity layers in the lower 
troposphere. From the degraded interference patterns the FSI algorithm underestimates 
the large bending angles produced by the vertical refractivity gradients leading to smaller 
refractivities and, ultimately, a negative bias with respect to ECMWF. 

We substantiate our interpretation of the difference between ZD and SD data with 
simulation results shown in Fig. 4. We note that the non-zero mean fractional refractivity 
error and standard deviation is caused by including L1/L2 carrier phase noise within the 
simulation. Above 20 km altitude the bending angles derived from the FSI are replaced by 
the values obtained from the geometric optics method, since the occurrence of multipath 
ray propagation in the stratosphere can be excluded. The simulations reproduce the over- 
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all differences between the observed ZD and SD refractivities (Fig. 3) as well as an (albeit 
smaller) bias between the two data sets at altitudes below 10 km. An additional SD 
simulation was conducted with the noise level of the reference link observations reduced 
to zero. Since the result (marked by crosses in Fig. 4) is in excellent agreement with the 
ZD profile throughout the full altitude range, the observed bias between the ZD and SD 
results in the lower troposphere is evidently caused by the reference link phase noise. 

4. Conclusion 

First radio occultation events observed by the GRACE-B satellite are successfully an- 
alyzed using zero differencing. In the upper troposphere and lower stratosphere the re- 
fractivities derived from zero and single differencing are in good agreement with the cor- 
responding ECMWF meteorological analyses. Zero differencing allows for a reduction in 
the amount of down-link data. More significantly, zero differencing reduces the noise level 
on the excess phase paths and thereby yields less-biased refractivities within regions of 
multipath signal propagation in the lower troposphere compared to single differencing. 
However, for quiet ionospheric conditions a modified single differencing method can be 
applied producing refractivities that are in excellent agreement with the zero differencing 
results. 
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Figure 1. Carrier phase noise extracted from the reference link observations as a function 
of voltage SNR V . The LI and L2 observations extracted from 109 reference link profiles are 
marked as circles and stars, respectively. The LI phase noise (full line) is calculated from Eqn. 6 
with T = 20 ms and B w = 20 Hz; the corresponding L2 noise (dashed line) is a heuristic 
parameterization since additional losses from semi-codeless or codeless L2 tracking are not taken 
into account by Eqn. 6. 
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Figure 2. Receiver clock drift derived from GRACE-B and CHAMP reference link data. Top: 
receiver clock drift observed in a CHAMP occupation event on 28 July 2004, 0:18 UTC at 69.4°N, 
12.5°W. Bottom: clock drift observed in a GRACE occultation event on 28 July 2004, 6:41 UTC 
at 66.6°S, 13.9°E (black). Note the change in vertical scale. 
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Figure 3. Left: Mean fractional refractivity deviation between GRACE observations and 
ECMWF analyses. The data set comprises 101 zero differencing (solid line) and 96 single differ- 
encing (dashed line) results. Thin lines indicate the 1-a standard deviations. Retrieval results 
from modified single differencing are marked as crosses. Right: number of data points retrieved 
as a function of altitude. 
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Figure 4. Same as Fig. 3, except that results from end-to-end simulations are shown. Left: 
the mean error for ZD and SD processing are plotted as thick solid and dashed lines, respectively. 
Thin lines indicate the l-a standard deviations. The result obtained from a SD simulation with 
noise-free reference link phase observations (crosses) is in almost perfect agreement with the ZD 
result. Right: number of data points retrieved as a function of altitude. 
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